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1. Introduction 


Induction heating is a physical process extensively used in the metallurgical industry for 
different applications involving metal melting. The main components of an induction heating 
system are an induction coil connected to a power-supply providing an alternating electric 
current and a conductive workpiece to be heated, placed inside the coil. The alternating 
current traversing the coil generates eddy currents in the workpiece and by means of ohmic 
losses the workpiece is heated (see Fig. 1). 

The design of an induction heating system mainly depends on its application. In this chapter 
we are interested in modelling the behavior of a coreless induction furnace like those used 
for melting and stirring. A simple sketch of this furnace is presented in Fig. 2(a). It consists 
of a helical copper coil and a workpiece formed by the crucible and the load within, which 
is the material to melt. The crucible is a cylindrical vessel made of a refractory material 
with higher temperature resistance than the substances it is designed to hold in. The coil, 
which is water-cooled to avoid overheating, is usually enclosed into a refractory material 
for safety reasons. Alternating current passing through the coil induces a rapidly oscillating 
magnetic field which generates eddy currents in the workpiece. These induced currents heat 
the load. When the load melts, the electromagnetic field produced by the coil interacts with 


Fig. 1. Induction heating furnace off (left) and on (right). Photographs courtesy of Mr. Victor 
Valcárcel, Instituto de Cerámica, Universidade de Santiago de Compostela. 
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Fig. 2. Sketch of the induction furnace (a) and radial section of inductors and workpiece (b). 


the electromagnetic field produced by the induced current. The resulted force causes a stirring 
effect helping to homogenize the melt composition and its temperature. 

It is well known that the high frequencies used in induction heating applications gives rise to 
a phenomenon called skin effect. This skin effect forces the alternating current to flow ina thin 
layer close to the surface of the workpiece, at an average depth called the skin depth. Thus, 
the ohmic losses are concentrated in the external part of the workpiece in such a way that the 
higher the frequency the thinner the skin depth. It is crucial to control the distribution of these 
ohmic losses, since they could cause very high temperatures in the crucible thus reducing 
its lifetime. Moreover, the frequency and intensity of the alternating current also affect the 
stirring of the molten bath. Since stirring mainly determines some of the properties of the 
final product, it is convenient to accurately know the influence of these parameters to achieve 
the desired result. 

Our goal is to understand the influence on the furnace performance of certain geometrical 
parameters such as the crucible thickness, its distance to the coil, the number of turns of the 
coil, or physical parameters such as the thermal and electrical conductivity of the refractory 
materials. In the last years, numerical simulation reveals as an important tool for this purpose, 
since it allows to introduce changes in the above parameters in silico, thus avoiding long and 
costly trial-error procedures in plant. 

The overall process is very complex and involves different physical phenomena: 
electromagnetic, thermal with phase change and hydrodynamic in the liquid region; all of 
them are coupled and it is essential to consider a suitable mathematical model to achieve a 
realistic simulation. 

Many papers have been published concerning the numerical simulation of induction 
heating devices, from some pioneering articles published in the early eighties (see, 
for instance, (Lavers, 1983) and references therein) to more recent works dealing with 
different coupled problems, such as the thermoelectrical problem appearing in induction 
heating (Chaboudez et al., 1997; Rappaz & Swierkosz, 1996), the magnetohydrodynamic 
problem related to induction stirring (Natarajan & El-Kaddah, 1999) and also a 
thermal-magneto-hydrodynamic problem (Henneberger & Obrecht, 1994; Katsumura et al., 
1996) but not fully coupled because material properties are not supposed to be dependent 
on temperature. Some other related works include mechanical effects in the workpiece 
(Bay et al., 2003; Hömberg, 2004). A more extensive bibliographic review can be found in 
(Lavers, 2008). 
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The aim of this chapter is to deal with the coupled thermo-magneto-hydrodynamic simulation 
of an induction heating furnace like the one described above. It is a survey of previous 
works by the authors (Bermúdez et al., 2007;b; 2009; Vazquez, 2008). Section 2 is devoted 
to describe the mathematical model to compute the electromagnetic field and the temperature 
in the furnace, and the velocity field in the molten region; the coupled mathematical model 
is a system of partial differential equations with several non-linearities. In Section 3a weak 
formulation of the electromagnetic model, defined in an axisymmetric setting, is derived. In 
Section 4 we present the time-discretization of the thermal and hydrodynamic models and 
the corresponding weak formulations. In Section 5 we describe the spatial discretization and 
the algorithms used to deal with the coupling formulation and the non-linearities. Finally, in 
Section 6, several numerical results are given for an industrial heating furnace designed for 
melting. 


2. Mathematical model for the induction furnace 


In this section we introduce the mathematical model to study the thermo-magneto-hydro 
dynamic behavior of an induction furnace. 


2.1 The electromagnetic model 

In order to model the electric current traversing the coil and the eddy currents induced in 
the workpiece, we introduce an electromagnetic model which is obtained from Maxwell 
equations. We state below a three-dimensional model defined in the whole space R? and 
deduce an axisymmetric formulation in Section 3. Since the full description from the 3D 
model to the axisymmetric one involves a lot of technical steps and it has been done in 
(Bermudez et al., 2009), we only give here a brief description and refer the reader to the quoted 
paper for further details. 

We begin by introducing some assumptions and notations related with the geometry that will 
be used in the sequel. Since the material enclosing the coil is not only a good refractory but 
also a good electrical insulator, the induction process in the workpiece is almost unaffected by 
the presence of this material. Thus, in the electromagnetic model, it will be treated as air. For 
simplicity, we do not consider the water refrigerating the coil. 

In order to state the problem in an axisymmetric setting, the induction coil has to be replaced 
by m rings having toroidal geometry. Let Op be the radial section of the workpiece and 
O4,O2,..-,Om be the radial sections of the turns of the coil, which are assumed to be simply 
connected (see Fig. 2(b)). Moreover, let us denote by Q the radial section of the set of 
conductors, i.e., inductors and conducting materials in the workpiece, given by, 


m 


A= UO, 
k=0 


and QS = R? \ ©. In particular, the refractory layer enclosing the coil is part of N°. 

Let A C RÊ be the bounded open set generated by the rotation about the z—axis of Q and 
A‘ the complementary set of A in RÌ. Notice that AS is an unbounded set corresponding to 
the air surrounding the whole device. Analogously, we denote by Ag, k = 0,...,m the subset 
of R generated by the rotation of Ox, k = 0,...,m, respectively, around the z—axis (see Fig. 
3). For simplicity, we will assume that Ao is simply connected and that A;, k = 1,...,m is a 
solid torus, in the sense that one single cutting surface Q% is enough for the set A, \ Ox to be 
simply—connected (see Fig. 3). 
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Fig. 3. Cylindrical coordinate system (a) and sketch of a toroidal turn A, (b). 


We denote by = the boundary of A and by I its intersection with the half-plane {(r,z) € R?;r > 
0}. We notice that & = Uj Xp, where X, denotes the boundary of Ay. Moreover, we assume 
that the boundary of Q is the union of T and Is, the latter being a subset of the symmetry axis 
(see Fig. 2(b)). 

The induction furnace works with alternating current so we can consider that all fields vary 
harmonically with time in the form: 


F (x,t) = Re [el F(x)], (1) 


where t is time, x € R is the spatial variable, i is the imaginary unit, F(x) is the complex 
amplitude of field F and w is the angular frequency, w = 27f, f being the frequency of the 
alternating current. 

For low and moderate frequencies, displacement currents can be neglected. Then, Maxwell’s 
equations can be reduced to the so-called eddy current model (see (Bossavit, 1998)): 


culH = J inR’, (2) 
iwB+curlE = 0 ink’, (3) 
divB = 0 inR, (4) 


where H, J, B and E are the complex amplitudes associated with the magnetic field, the current 
density, the magnetic induction and the electric field, respectively. 
System (2)-(4) needs to be completed with the following constitutive relations: 


B = uH inR’, (5) 
cE inA 

ne á 6 

J f in AS, (6) 


where y is the magnetic permeability -which is supposed to be independent of H- and ø is 
the electric conductivity. Both are assumed to be bounded from below by a positive constant 
and dependent on both the spatial variable and the temperature T, i.e., y = u(x, T) and g = 
a(x,T). Nevertheless, for the sake of simplicity, we drop this dependency in the notation until 
forthcoming sections. 
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Notice that we have neglected the fluid motion in the Ohm’s law (6), because a dimensionless 
analysis shows that it is negligible in comparison with the other terms for this kind of 
applications. 
Since the previous equations hold in R, we also require for the fields a certain behavior at 
infinity. More specifically, we have, 
E(x) = O(|x|~*), uniformly for |x| — œ, (7) 
H(x) = O(|x|7!), uniformly for |x| — oo. (8) 
Model (2)-(8) must be completed with some source data related to the energizing device. 


In particular, we would like to impose the current intensities I = (l1,...,Im) crossing each 
transversal section of the inductor, i.e., 


| J-t=h, k= tpn, (9) 
Ok 
where v denotes a unit normal vector to the section Ox. 


Remark 21 Imposing these constraints is not trivial at all and requires to relax in some sense equation 
(3). We refer the reader to Remarks 2.1 and 2.2 of (Bermúdez et al., 2009). To this respect, we also 
cite the paper (Alonso et al., 2008) where the authors give a systematic analysis to solve eddy current 
problems driven by voltage or current intensity in the harmonic regime and in bounded domains. 


To solve the eddy current problem in an axisymmetric setting, we are going to propose a 
formulation based on the magnetic vector potential. To do that, we also need to introduce a 
suitable scalar potential. Firstly, equation (4) allows us to affirm that there exists a magnetic 
vector potential A defined in R? such that, 


B = curl A, (10) 
and from equation (3), we obtain 
iwcurlA+curlE=0 inA. 


Taking into account the form of the kernel of the curl operator in each connected component 
of the conductor, we can say that (see, for instance (Amrouche et al., 1998)) 


iwA+E=-—v ind, (11) 
where 
v= gradU, 


and U is a scalar potential having a constant jump through each Ox, for k = 1,---,m. We have 
denoted by grad the gradient operator in the space H! (A \ U"_,Q;). As we will show below, 
this representation allows us to impose the sources in the closed circuits A, (see also Section 
5.2 of (Hiptmair & Sterz, 2005)). 

For k = 1,...,m, let us denote by yx the solution in H! (Ag \ Ox), unique up to a constant, of 
the following weak problem: 


f ogradn,-gradé=0, Vē € H! (AŅ), (12) 
Ar \ Ok 


[mlo = 1, (13) 
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where [7;]o, denotes the jump of 7; through Ox along v. 
By using functions g, k = 1,---,m, the scalar potential U can be written as (see again 
(Amrouche et al., 1998)), 


m 
U=0+ )° Vin, (14) 
k=1 


with ® € H! (A) and yk =1,...,m being some complex numbers. From the definition of 7, 
we deduce that V; is the constant jump of ü through each surface Ox, k = 1,---,m. From 
a physical point of view, these complex numbers V; can be interpreted as voltage drops 
(Hiptmair & Sterz, 2005). 


Then, taking into account that H = u~! 


curl A and equation (2) we obtain 


m 
iwoA + curls curl A) = —ov = —o(grad® + } ` Vi. gradi). (15) 
k=1 

Notice, however, that the vector potential A is not unique because it can be altered by any 
gradient. Thus, in order to get uniqueness we need to impose some gauge conditions. In the 
conductor region A, we set div(7A) = 0 and the boundary condition 7A - n = 0 on £, where 
n is a unit normal vector to & outward from A. 
From these gauge conditions, we can easily deduce that grad ® = 0. Consequently, 


m 
v= } grad ny 
= 


and equation (15) can be written as 


m 
iwoA + curi(7 curl A) = —ov = -0 L Vegfadyę ind. (16) 
k=1 


Notice that, in particular, the electric field in A is given by 


m 
E=-—iwA-— ) V gradn,. (17) 
k=1 


On the other hand, in the air region, we impose the gauge conditions 
divA=0inA° and | A-n=0, j=0,...,m. (18) 
Xj 


In Section 3we will obtain a weak formulation from equations (16), (18) and the imposition of 
the current intensities (9) in the axisymmetric setting. 


2.2 The thermal model 

The electromagnetic model must be coupled with the heat equation to study the thermal 
effects of the electromagnetic fields in the workpiece. We will describe the equations of the 
model in an axisymmetric setting, paying attention to the terms coupling the thermal problem 
with the electromagnetic one. 
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Symmetry axis 


Fig. 4. Computational domain for the thermal problem 


The computational domain for the thermal model is the radial section Qr of the whole furnace, 
that is, 
Ox = Qo U Q1 U- U Qm U Ont 


where ©,,,41 denotes the radial section of the dielectric parts of the furnace (see Fig. 4). We 
notice that we now take into account that the coil is water-cooled and replace each turn by a 
hollow torus. An appropriated boundary condition will be considered on the inside boundary. 


Remark 22 For simplicity, in Section 2.1 the coil was replaced by rings with toroidal geometry, i.e., 
by solid tori. Since we are now considering a refrigeration tube along the coil, it will be replaced by 
hollow tori. We refer the reader to Remark 4.4 of (Vazquez, 2008) to see that, under axisymmetric 
assumptions, the electromagnetic model does not change in this topology. 


Since the metal is introduced in solid state and then melted, we use the transient heat transfer 
equation with change of phase, written in terms of the enthalpy. Furthermore, since the 
molten metal is subject to electromagnetic and buoyancy forces, we also need to consider 
convective heat transfer. Let us suppose that we already know the velocity field u which is 
null in the solid part of the workpiece, and the current density J. Then, the equation for energy 
conservation is 


2 
é — div(k(x,T) grad T) = ug in Or, (19) 


where T is the temperature, k is the temperature-dependent thermal conductivity and é is the 
material derivative of the enthalpy density. It is given by 


__ de 
Ot 


é +u- grade. 


The source term on the right-hand side of (19) represents the heat released due to the Joule 
effect; it is obtained by solving the electromagnetic model introduced in the previous section. 


Remark 23 Notice that the time scale for the variation of the electromagnetic field is much smaller 
than the one for the variation of the temperature. Indeed, the physical parameters used in a typical 
industrial situation give a time scale for temperature of the order of 60 seconds while for the 
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electromagnetic problem it is of the order of 10~> seconds. Thus, we may compute the electromagnetic 
field in the frequency domain using the eddy current model and then the heat source due to Joule effect 
is determined by taking the mean value on a cycle, namely, 
w 2n/w 
Za ha I (x,t) E(x, t)dt, (20) 
J and E being the current density and the electric field, respectively, written in the form of (1). It is 
not difficult to prove that expression (20) coincides with the right-hand side of (19). 


In general, the metal in the crucible undergoes a phase change during the heating process. For 
this reason, the enthalpy in (19) must take into account the latent heat involved in the phase 
change. This can be done by introducing an enthalpy function similar to that used for Stefan 
problems (see (Elliot & Ockendon, 1985)). More precisely, the enthalpy density is expressed 
as a multi-valued function of temperature by 


e(x,t) E H(x,T(x,£)), (21) 
with 
¥(x,T) T< T(x), 
H(x,T)= 4 [¥(x,Ts),¥(x,Ts) + L o(x,Ts)] T = T,(x), (22) 
¥(x,T) +L p(x,Ts) T > Ts(x), 
and 
T 
YT) = f pge?) ae, (23) 


p being the density, c the specific heat and L the latent heat per unit mass, i.e., the heat per 
unit mass necessary to achieve the change of state at temperature Ts. 


Remark 24 It is to be noted that H is a multivalued function rather than discontinuous. Indeed, 
an element of volume at the solidification temperature may have any enthalpy value in the interval 
[F(x Ts), P(x, Ts) + p(x, Ts) L]. 


2.2.1 Thermal boundary conditions 

Equation (19) must be completed with suitable initial and boundary conditions. Let Tr be the 
boundary of Or; we denote by I's its intersection with the symmetry axis, by Ig the part of the 
boundary in contact with air and by Tc the internal boundary of the coil which is in contact 
with the cooling water (see Fig. 4). On Tc we consider a convection condition 


oT 
on 


a being the coefficient of convective heat transfer and Twy the temperature of the cooling water. 
On boundary Ip we impose the radiation-convection condition 


k(x, T) = a( Tw = T), (24) 


k(x T) or =a(T; —T)+7(T —T*), (25) 


an 
where T; and T, are the external convective and radiative temperature, respectively, and y is 
the product of the emissivity by the Stefan-Boltzmann constant (5.669e-8 W/ m°K*). Finally, 
on I's we set the symmetry condition 


k(x,T)— =0. (26) 
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Molten metal 


Q, 


$7] 


Fig. 5. Computational domain Q} (t) for the hydrodynamic problem. 


2.3 The hydrodynamic model: Boussinesq approximation 

As mentioned before, in order to achieve a realistic simulation of the overall process occurring 
in the furnace, convective heat transfer must be taken into account. The hydrodynamic 
domain is the molten region of the metal, which varies as the metal melts or solidifies, 
making our hydrodynamic domain time-dependent. This molten metal is subjected to both 
electromagnetic and buoyancy forces, the latter due to the variation of mass density with 
temperature. 

Since the range of temperatures in the molten region is not very large, we use the Boussinesq 
approximation to model the fluid motion. Roughly speaking, this approximation consists on 
assuming that the mass density is constant in the inertial term, and that it depends linearly 
on the temperature in the right-hand side. Denoting by 0,(t) the radial section of the molten 
metal, and by s(t), T4(t) and Tn (t) the different parts of the boundary at time t (see Fig. 5), 
the fluid motion is modeled by the following equations 


ou e 
po (St + (gradu)u ) — div2yoD(u)) + grad p = pbt- Toetie 27) 
divu = 0, (28) 


where u is the velocity field, p is the (modified) pressure, and D(u) is the strain rate tensor, 
namely D(u) = (gradu + gradu’)/2. Constants pọ, yo and ßo denote the mass density, 
the dynamic viscosity and the thermal expansion coefficient, respectively, at the reference 
temperature To. The first term in the right-hand side takes into account the buoyancy forces (g 
denotes the gravity acceleration), whereas f; represents the Lorentz force which is computed 


i w [27/w 
fi = a Jp J (x,t) x B(x,t) dt, (29) 


where B is the magnetic induction field written in the form (1). 
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Equations (27)-(28) are completed with the following boundary and initial conditions: 


u =0 onT4(t), (30) 
Sn =0 onTI;(t), (31) 
Sn =0 onT;(t), (32) 

u =0 inQ,(0), (33) 


where S denotes the Cauchy stress tensor, S = — pI + 279D(u) and I is the identity tensor. 
When using the Boussinesq approximation, the thermal model is modified in such a way that 
the material properties in the molten region are considered at the reference temperature, that 
is to say, 


PT ul gradat) -divttygraa y -UÈ (34) 
Poco > +4 gra iv(ko gra Noe 
where co and kp represent the specific heat and the thermal conductivity at the reference 
temperature, respectively. We remark that this approximation is only used in the molten 
region of the metal. In the rest of the domain the heat equation remains non-linear. 


2.3.1 An algebraic turbulence model: Smagorinsky’s model 
To take into account the possibility of turbulent flows, it is necessary to introduce a turbulence 
model. The basic methodology consists in modifying equations (27), (28) and (34), in the form 


ot . ` z R 
Po (= + (grad) 4) — div(2i7e¢,D(a)) + gradp = —poßo(T — To)g +f, (35) 
diva = 0 (36) 
dð r ; A J|? 
poco | 3, +û-gradT ) — div(keffgrad T) = oer (37) 


where he ff and ke ff denote the effective viscosity and thermal conductivity, respectively; 
namely, 17, ff = 10 + Nt and k, ff= ko + kt, 4¢ being the turbulent viscosity and k; the turbulent 
thermal conductivity. 

Notice that fields u, p and T in (27), (28) and (34) have been replaced by their corresponding 
filtered fields û, p and T (see (Wilcox, 2008)). Hereafter, and for the sake of simplicity, we will 
eliminate the symbol * when referring to these variables. 

For a rigorous derivation of turbulence models we refer the reader to, for instance, 
(Mohammadi & Pironneau, 1994; Wilcox, 2008). The turbulence models differ essentially on 
the way the turbulent viscosity 7; and the turbulent conductivity k; are computed. An efficient 
and easy to implement model is the one proposed by Smagorinsky, which has been considered 
in the present work. It belongs to the family of Large Eddy Simulation (LES) models. In this 
case, 7; and k; are given by 


nt =poCh2|D(u)|, C 0.01, k= oot, (38) 
t 
where h(x) is the mesh size of the numerical method around point x, and Pr; is the turbulent 
Prandtl number, which is taken equal to 0.9. 


It is worth noting that more accurate turbulent models, such as the well-known k-e model, 
can also be used. The drawback is that the computational effort to apply those models is 
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much higher, since they require the solution of some additional partial differential equations 
in order to compute the turbulent viscosity y+. However, if an accurate simulation of the fluid 
motion is necessary -as for instance in electromagnetic stirring- the use of one of these models 
might become unavoidable. 


3. Weak formulation of the electromagnetic problem in an axisymmetric setting 


In this section we present the weak formulation of the electromagnetic problem. We start by 
introducing some functional spaces and sets. Let X be the Beppo-Levi space (see (Nédélec, 
2001), Section 2.5.4) 


G(x) 


 /1+ |x? 


X4{G € L?(R?), curlG € L?(R3)}, 


and its subset 


Y ={GEX:divG =0in ds, | G-n=0,j=0,...,m}. 
xj 


We recall that we are interested in finding a solution of the eddy current problem satisfying 
the intensity conditions (9). To attain this goal, we start by writing these constraints in a weak 
sense. 

Since current density J = cE satisfies div J = 0 in A and J -n = 0 on È, we have for k = 1,...,m, 


Jeeradm = -f  divimt+ fo Joam f m 69) 
i grad 1 ae TTT lt 
where vz is the unit vector normal to the radial section OQ. Hence, we can impose the current 
intensities as follows: 


m m 


Laf cE: grad ýk = L We YW = (W1,..., Wm) € C”. 
kzi YAA% k=1 


Then, taking into account (17), we obtain the following weak form of constraint (9) which is 
well defined for any vector function A € Y and vector W € C”: 


- wf iwogiadn,- A- 
Ł TN i 


On the other hand, multiplying equation (16) by the complex conjugate of a test function G, 
denoted by G, integrating in IR° and using a Green’s formula we can easily obtain, 


m m 


Lw i oV;,| grad ny? = > We. (40) 
kat An \O% k=1 


2 1 2 ul = 
iw vA-G+ f —curlA-curlG + vf ogrady,-G, VGEY. 41 
i ke er” L OE y (41) 
Thus, we are led to solve the following mixed problem: 


Problem MPI.- Given I = (lh,...,lm) € C”, find A € Y and V € C", such that 


m 


wf cA-G+ f TE e Luj ogrady,-G=0,VGE V. 
JIR3 R? H k Ak\ Ox 


=] 
m 1 2 2 1, Z 
wf opiady, A+ — wy | oV;| giad]? = —— ) Wey VW EC”. 
L kag Pa wa STA k| grad y| uS kWk 
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Notice that the vector field V of voltage drops, can be interpreted as a Lagrange multiplier 
introduced to impose the current intensities in a weak sense. The mathematical analysis of 
this mixed formulation can be found in (Bermúdez et al., 2009), where the simpler formulation 
with the voltage drops as data is also studied. 


3.1 An axisymmetric BEM/FEM formulation of problem MPI 

Notice that the above formulation is written in the whole space RÌ; in (Bermúdez et al., 
2007) we have approximated the problem in a bounded domain by defining approximated 
boundary conditions far from the conducting region and using a finite element technique. 
However, in this work, we consider a hybrid boundary element/finite element method (in 
the sequel BEM/FEM) to solve the problem in the whole space. To attain this goal, we are 
going to write the problem MPI in another form involving only the values of the magnetic 
vector potential A in A and on its boundary =. Notice first that the field pw} curl A, which is 


the intensity of the magnetic field, belongs to 1, and then its tangential trace wt curlA xn 
is continuous across &. Besides 
1 
curl (+ curl a) =curlH=0 inAS, (42) 
0 


where [ig denotes the vacuum magnetic permeability. Then, by using a Green’s formula in A‘, 
we have, 


Sa -curlG =f = curl A -curlG +f 1 curl A - curlG 
JR 4 AH Ac Ho 
1 E 1 - 
= f ~curlA-curlG — | —curlA xn-G,VGEy. 
JAH x Ho 
Thus, the first equation of problem MPI can be formally written as: 


iw [oa G+ f Scart cue — | cane ae 
JA JAH x Ho 


m 


Evf ograd: Č =0, YG € Y. (43) 
k= IANO 


We notice that the value of ug leurl A x n on E can be determined by solving an exterior 
problem in A‘. Since we are interested in the numerical solution of the problem in an 
axisymmetric domain, we directly transform this term in the axisymmetric case. 

We consider a cylindrical coordinate system (r,0,z) with the z—axis coinciding with the 
symmetry axis of the device, (see Figure 3). Hereafter we denote by er, eg and ez the local 
unit vectors in the corresponding coordinate directions. Now, cylindrical symmetry leads us 
to consider that no field depends on the angular variable 6. We further assume that the current 
density field has non-zero component only in the tangential direction eg, namely 


J= Je(r,z) ep. 


We remark that, due to the assumed conditions on J, (3), (6) and (10), only the @-component of 
the magnetic vector potential, hereafter denoted by Ag, does not vanish, i.e., 


A = Ao(r,z)eg. (44) 
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Thus A automatically satisfies (18). Let G = (r,z)eg be a test function and n = n; e; + nz ez. 
By taking into account the expression for curl in cylindrical coordinates and the fact that 


1 
grady = Iy” in Ay, k=1,...,m, (45) 
the axisymmetric version of the problem formally writes as follows: 


Given I = (Iy,...,lm) € C”, find Ag and V € C", satisfying, 


iw f Ag- prdrdz 4 a i ae | ROA Ade 


ne or ðr Qu Oz dz 
1 SAG); B 7 
Tbh oa mÈ Me fl opdrdz =0, (46) 
1 # a 1 m W 7 
32 Ag drdz)W, — drdz)W, IW, 47 
ee MONG 9 drdz) eae fie! rdz)W, = Sti bh We, (47) 


for all test function p and W € C”. 


On the other hand, the term feug 1a(rAg)/anpdy can be transformed by using the 
single-double layer potentials. We refer the reader to (Bermúdez et al., 2007b) for further 
details concerning this transformation and introduce the same notation of that paper, namely 


Ag = rAg, 
i dA, dA, 
A (rz) = g + a z 


We are led to solve the following weak problem 


Problem WEPI.- Given I= (Å,..., Im) € C”, find Ap :Q— C, V € C” and X :T — C such that 


A CO gr al $ =) 1 r=; T Z CO -, = 
wf T Agy drdz + | -p erad Ay -grady drdz — fA pdyt+ Zn Enja? drdz = 0, 


1 m = C med 1 m _ dd ji m 
— Yw ZA W -k drdaz = -— Y I 
Ir Oe ef? garde + ao ‘foo ma Qrtiw dL w 


k=1 


| r 1 los rae 
[Lager ihata [Morten 


for all test functions y, Ç and W € C™. In the previous equations, G and Gn denote 
the fundamental solution of Laplace’s equation and its normal derivative, respectively, in 
cylindrical coordinates (see again (Bermúdez et al., 2007b)). 


4 Time discretization and weak formulation of the thermal and hydrodynamic 
problems 


In this section we introduce a time discretization and a weak formulation of the thermal and 
hydrodynamic models. Again, we exploit the cylindrical symmetry and we assume that u 
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does not depend on @ and it has zero component in the tangential direction eg. In order to 
simplify the notation, in what follows we shall drop index t for O; (t). 

To obtain a suitable discretization of the material time derivative in equations (19) and (27) we 
shall use the characteristics method (see, for instance, (Pironneau, 1982)). 

Given a velocity field u we define the characteristic curve passing through point x at time t as 
the solution of the following Cauchy problem: 


d 
qk tT) = (Xx, 7),7), (48) 
X(x,t;t) =x. 


Thus X(x,t;T) is the trajectory of the material point being at position x at time t. In terms of X, 
the material time derivative of e is defined by 


0t) = E EXT), T) (49) 


Let us consider a time interval [0,t/] and a discretization time step At = tf/N to obtain a 
uniform partition II = {t” = nAt,0 < n < N}. Let e” and u” be the approximations of e and u 
at time t”, respectively. We approximate the material time derivative of e at time pnt by 


6(x, pert) ~ en (x) =el (x"(x)) 


— ; (50) 


where x” (x) = X” (x,t"+1;t") is obtained as the solution of the following Cauchy problem 


d yn (xa 5c) = u(x" (x, Vl rhe), 661) 
X"(x, ott, ntl) =x, 

at time t”. Notice that, since u = 0 in the solid region, the solution of this Cauchy problem is 
X”(x,t" t1; T) = x for any T. Hence equation (50) in the solid part is equivalent to a standard 
time discretization without using the method of characteristics. 

Analogous to (50), we consider in (35) the above two-point discretization for the material time 
derivative of the velocity. Thus, taking into account the cylindrical symmetry, multiplying 
equations (35) and (36) by suitable test functions, and integrating in the liquid domain 
O; we obtain, after using the Green’s formula, the following weak formulation of the 
semi-discretized hydrodynamic problem (the : representing the scalar product of two tensors) 


Problem WHP.- For each n = 0,1,...,N — 1, find functions u"*? and p"*+ such that u"*! = 0 on 
Tq and furthermore 


1 


me pou”! - wrdrdz + I, Nef f (gradu"*? : gradw) rdrdz + 


L Nef f ((gradu"*1)! : grad w) rdrdz — h p"* div wrdrdz = 


1 
— T” — Ty)g - wrdrdz +f £°+1. wrdrdz + zl "1 ox") . wrdrdz, 
f eobot 0)8 a Mi am o x") 


f divu” t! grdrdz = 0, 
O 
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for all test functions w null on T4, and q. 


On the other hand, assuming cylindrical symmetry and the fact that the temperature does not 
depend on the angular coordinate 0, we can write equation (19) in cylindrical coordinates. 
Then, applying the time discretization (50), multiplying by a suitable test function and using 
a Green’s formula we obtain the following weak formulation of the semi-discretized thermal 
problem: 

Problem WTP.- For each n = 0,1,...,N — 1, find a function T"+1 such that 


Í, Z e+! Zrdrdz + T kegp(t,z,T"*1) grad T"*1 . grad Zrdrdz 
JOr : 


L l a(To SB Zral + | («(Te — T!) + of Te — (T"+1)4))Zrdr 
JTe JTR 


1 n n L 1 +1)2 
+ l aa ox") Zrdrdz 4 Oo, Jone, TA) go | Zrdrdz, 


for all test function Z. 


Remark 41 Note that, from equations (6), (17) and the expression for functions yg, we can infer that 


Jo = —iwoAg in Qo, 
V, 
g = -iwr Ag — = m Ory k=1,...,m, 
27r £ 


Jo =0 in Om+1- 


Hence, to determine the heat source we need to compute an approximation of field Ag at time "+t. 
This is determined as the solution of the weak formulation WEPI with the physical parameters u and 
o being evaluated at temperature T’*1. 


5. Space discretization and iterative algorithm 


Problem WTP has been spatially discretized by a piecewise linear finite element method 
defined in a triangular mesh of the domain Qr. For the spatial discretization of problem WEPI 
we have used a BEM/FEM method (see (Bermúdez et al., 2007b) for further details). Finally, 
problem WHP has been spatially discretized by the finite element couple P4-bubble / P4, which 
is known to satisfy the inf-sup condition (Brezzi, 1991); this last problem only takes place in 
the time-dependent liquid domain Q; (t), which must be computed at each time step. 

It is important to notice that, at each time step, there are several terms coupling the three 
problems. However, since we are neglecting the velocity in Ohm’s law, and the method of 
characteristics is used with the velocity at the previous time step, the hydrodynamic problem 
can be solved uncoupled from the two others. Nevertheless, the coupling between the thermal 
and the electromagnetic models cannot be avoided: the heat source in the thermal equation 
is the Joule effect and the solution of the electromagnetic problem depends on the electrical 
conductivity and the magnetic permeability, which may vary with temperature. 

Moreover, the thermal problem WTP contains several nonlinearities: 


— The thermal conductivity k depends on temperature. 


— The external convective temperature Tw also varies with the heat flux of the inner 
boundaries of the coil. 
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— The enthalpy e depends on temperature and it is a multivalued function. 


— The radiation boundary condition depends on T$, 


To treat the nonlinear terms we are going to introduce several iterative algorithms. The 
dependence of the thermal conductivity k can be easily treated, just by taking in the heat 
equation the thermal conductivity evaluated at the temperature of the previous time step (or 
of the previous iteration of an outer loop). This can be done because the thermal conductivity 
is a smooth function. Nevertheless, as we will show in the next section, the other nonlinearities 
need other more sophisticated iterative algorithms to be introduced. 


5.1 Iterative algorithm for the temperature of cooling water 

As we said before, the induction coil of the furnace is water-cooled to avoid overheating; this 
fact is modeled by means of boundary condition (24), where the temperature of the cooling 
water, Tw, depends on the solution of our problem which introduces a nonlinearity in the 
equations. 

To deal with this nonlinearity we will seek the convergence of the heat flux from the coil to 
the cooling water. From the solution of the thermal problem this heat flux is computed as 


aT 
H=2n f kgr, (52) 


where the integral is multiplied by 271 because Tc is only a radial section of the boundary. 
Using (52) the outlet water temperature can be computed as 


penie; 63) 
PwlwQ 

with pw and cw the density and specific heat of water, respectively, T; the inlet temperature of 
the cooling water, and Q the water flow rate. 
To solve the problem, we assume that water temperature Tw is constant along the coil which 
is actually a reasonable assumption, since the difference between the inlet and the outlet 
temperature is seldom higher than 10 °C. 
The temperature of the cooling water is then computed using an iterative algorithm. Let us 
suppose Tw,j is known. Then, at iteration j + 1, we compute T;,; as the solution of thermal 
problem WTP with Tw = Tw, j. Then, we compute the heat flux H from equation (52) and set 


H 
20wlwQ i 


ie., Tyj41 is the mean value of the given inlet temperature and the computed outlet 
temperature. 

It is worth noting that an implicit algorithm is needed to avoid instabilities. The iterations of 
the implicit method can be merged with the iterations of the thermoelectrical coupling leading 
to a low computational cost. 


Tw,j+1 = T; Bhs 


5.2 Iterative algorithms for the enthalpy and the radiation boundary condition 

In order to solve the non-linearities due to the multi-valued character of the enthalpy and 
to the exterior radiation boundary condition, we propose a fixed point algorithm which is 
described in detail in references (Bermúdez et al., 2003; 2007; Vazquez, 2008). It is summarized 
here for the reader’s convenience. 
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At time step (n + 1) and for a positive 6 we introduce the new function 


ie = ett _ pre", 


Using (21) and the fact that H is a maximal monotone operator, we get 


g** (7,2) = HE (1,2), T**1(r,z) +Ag™*1(r,z)) for O0<A< ie (54) 


uh being the Yosida approximation of the operator HP = H — BT (Bermtidez-Moreno, 1994). 
The same idea is introduced to deal with the nonlinearity associated to the fourth power of 
the boundary temperature in (25): we consider the maximal monotone operator G(T) = |T|T°, 
which coincides with T4 for T > 0, and we define the new function 


gttl = G(T"*?) } cT'tt, 
Then using G5, the Yosida approximation of operator G% = G — xT , we get 


se gy ONT Gg) + ôs”+1(r,z)) for 0< ô< a (55) 


Now, in order to solve problem WTP, the idea is to replace e”+1 by g”+1 + BT"*! and G(T"*") 
by s"*1 + x T”+1, Finally, to determine T”+1, q+! and s"+! we introduce an iterative process 
using (54) and (55). 

We notice that the performance of the proposed algorithm is known to depend strongly on 
the choice of parameters £, A, ô and x. In (Vázquez, 2008) an automatic procedure is proposed 
to compute these parameters as functions of (r,z), which accelerates the convergence of the 
method. 


5.3 Iterative algorithm for the whole problem 

Now we present the iterative algorithm to solve the three coupled models, along with 
the nonlinearities. Basically, it consists of three nested loops: the first one for the time 
discretization, the second one for the thermoelectrical coupling, and the third one for 
the Bermudez-Moreno algorithms for the enthalpy and the radiation boundary condition 
presented above. As we have said before, the hydrodynamic problem is solved at each time 
step, uncoupled from the two others. Moreover, the nonlinearities of the thermal conductivity 
k and the cooling water temperature Ty are also treated by iterative algorithms and their 
corresponding loops are in fact merged with that of the thermoelectrical coupling. For the 
sake of simplicity we present in Fig. 6 a sketch of the algorithm with the three nested loops. 


5.3.1 The algorithm 

Let us suppose that the initial temperature T? and velocity u? are known. From T° we 
determine the initial enthalpy e° and set the temperature of cooling water T°, = T;, being 
T; the inlet temperature. Then, at time step n, with n = 1,...,N we compute Aj, T” and u” by 
doing the following steps: 


1. If u”! £0, compute x"~!(x), the solution of (51). 


2. Calculate the turbulent viscosity 4# = poCh?|D(u"—!)| and the turbulent thermal 
conductivity ký = coy}! / Prt. 
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YES 


Fig. 6. Scheme of the numerical algorithm considering all the nested loops. 
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3. Set T = ToL eg =e" 1 and Too = Tw 1. Then compute Ay, T”, e” and Th as the limit of 
Aj j T! and T” j obtained from the following iterative procedure: 
1. For j > 1 let us suppose that Tia and Toj- , are known. Then set Aj = = Ap /r, with Ae 
A’ and V € C” the solution of WEPI by taking o = o(r,z „T? ,)and u = p(r,z,T! , 
y 8 j H=H j 
(b) Set Jf), = —iwo(r,z,Tj' 1) Ag; — Vk/ (27r) in Ok, k = 0,...,m (with Vo = 0). 
(c) Determine the optimal parameters x«'’(r,z) and f"(r,z) with the method described in 
P P j j 
(Vázquez, 2008). 
(d) Set Tjo = Ti and also 
Ijo — -p Tja 
So = te = “ett, 
Then compute T, q and s7 as the limit of the following iterative procedure: 


, are known. Then compute Tik as the 


1. For k > 1 let us suppose that qj x1 and sh k- 


solution of the linear system 


A Lp " Zrdrdz + f kT 1) +K?) grad T7 k gradZrdrdz + | a Zr al 


+f (a +y?) jel ra =f aT j- 1Zrdr ei (aTe + yTe — Sj 4—-1)Zr aT 
Tr í Ir í 


iene 4 1 2 : 
+ eae ox” -= qi) Zrdrdz + Ja EZT M Zrdrdz VZ test function. 


ii. Update multipliers Iik and Sik by means of the formulas 


Vik = Hi ( Tik t Agia i), 
Sik = G(T + 6874-1). 
(e) Update the value of the enthalpy and the value of the cooling water temperature by 
computing 
ni a npn 
TT 
oT; i H 
H=2 k—rdf, d T3 , = T; +=. 
T Tce on r a Or] it 20wlwQ 


4. Determine the hydrodynamic domain from the value of enthalpy e”. 
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5. Find u” and p” solution of 
1 
KI a pou" - wrdrdz + J, Cr +ý) (gradu” : grad w) rdrdz 
+ (yo +47) ((gradu")! : gradw) rdrdz — i p” divwrdrdz = 
l l 
-f poßbo(T” — To)g : wrdrdz + k £,(Aj) : wrdrdz + 
l 1 
af polu”! o x"7!)-wrdrdz, V test function w, 


f divu"q=0, V test function q. 
Q; 


6. Numerical simulation of an industrial furnace 


In this section we present some numerical results obtained in the simulation of an 
industrial furnace used for melting and stirring. The results have been performed with 
the computer Fortran code THESIF (http://www.usc.es/~thesif/) which implements the 
algorithm described above. 


6.1 Description of the furnace 

We briefly describe the geometry and working conditions of the furnace and refer the reader 
to Chapter 4 of (Vázquez, 2008) for further details. 

The inductor of the furnace is a copper helical coil with 12 turns which contains a pipe 
inside carrying cool water for refrigeration. Inside the coil a crucible is placed, containing 
the metal to be melted. The crucible is surrounded by an alumina layer to avoid heat losses. 
The reason to take alumina for this layer is that it is not only a good refractory but also 
a good electrical insulator. Thus the induction process in the crucible and in the metal is 
almost unaffected by the presence of alumina. For safety reasons, the induction coil is also 
embedded in the alumina layer. Above this alumina layer there is a layer of another refractory 
material, called Plibrico. Moreover, for safety reasons there is a metal sheet surrounding the 
furnace, which prevents from high magnetic fields outside the furnace. The furnace rests 
on a base of concrete, which will be also considered in the computational domain. For the 
same purpose of thermal insulation, there is also placed a lid over the load, but it will not be 
considered it in the numerical simulation. Otherwise, we should solve an internal radiation 
problem, instead of imposing the convection-radiation boundary condition given in (25). 
In (Bermúdez et al., 2006; 2011) an improved thermal model including an internal radiation 
condition is introduced. 

The computational domain for the electromagnetic problem consists of the furnace (load, 
crucible and coil) as it is shown in Fig. 7. The computational domain for the thermal problem 
is composed by the metal, the crucible, the refractory layers and the coil. The computational 
domain for the hydrodynamic problem is the molten metal and it is computed at each time 
step. 

Concerning the values of the physical properties of the different materials, we show only the 
electrical conductivity of the crucible and metal to be melted because it plays a major role 
in the results. We notice that the electrical conductivity depends on temperature, being this 
dependency really important in the second case, because the load is assumed to be an insulator 
when solid and a good conductor when it melts, as can be seen in Fig. 8. 
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Fig. 7. Distribution of materials in the radial section of the furnace. 


In the electromagnetic model that we have presented two working parameters have to be 
chosen: the frequency of the alternating current and the corresponding current intensity. 
However, when trying to simulate the real furnace, the given data is the power, and the 
current intensity is adjusted to obtain that power in the furnace. Moreover, since the electrical 
conductivity of the materials varies with temperature, the relationship between power and 
intensity also changes, and the intensity is dynamically adjusted during the process. To deal 
with this difficulty, the algorithm was slightly modified (see (Bermúdez et al., 2011) for further 
details) to provide the power as the known data and then to compute the intensity to attain 


the given power. Thus, for the numerical simulation, we have used the total power of the 
system as data. 


6.2 Numerical results 


We have performed two numerical simulations of the furnace with the same value of the 
power and two different values of the frequency, 500 Hz and 2650 Hz, to see how the 
frequency affects the heating and stirring of the metal. 

In Fig. 9 we represent the temperature in the furnace for each simulation. As it can be seen 
the temperatures obtained in the furnace are very similar, but a little higher in the case of the 
highest frequency. We notice the strong influence of the refrigeration tubes in the temperature: 
the temperature in the copper coil and the surrounding refractory is about 50 °C, causing a 
very large temperature gradient within the refractory layer. In Fig. 10 we show a detail of 
the temperatures in the crucible and in the load. As it can be seen, higher temperatures are 


5 5 
25%10 4g X 0 


N 


Electrical conductivity (27' m™') 
Electrical conductivity (7! m~) 


t 0 
0 500 1000 1500 2000 0 500 1000 1500 2000 
Temperature (°C) Temperature (°C) 


Fig. 8. Electrical conductivity for the crucible (left) and the load (right). 
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Fig. 9. Temperature after three hours for 500 Hz (left) and 2650 Hz (right). 


reached when working with high frequency. This fact is explained due to the distribution of 
ohmic losses, as we will explain below. 

In Fig. 11 we show the Joule effect in the load and the crucible after five minutes, when the 
metal is still solid. In Fig. 12, the same field is represented after three hours, when the metal 
has been melted. Comparing the results for different frequencies, we can see that the higher 
the frequency the higher the maximum values of the Joule effect, but due to the skin effect 
they are concentrated on the external wall of the crucible. Decreasing the frequency allows 
a better power distribution, at the cost of using higher intensities, thus causing larger power 
losses in the coil (see Fig. 14). 

Moreover, comparing the results for solid and molten load, we see how the high conductivity 
of the molten metal affects the performance of the furnace. At low frequency the ohmic 
losses in the load become higher when the material melts, thus heating the load directly 
and reducing the crucible temperature (see Fig. 14 for the ohmic losses and Fig. 10 for the 
crucible temperature). Moreover, the presence of molten material increases the skin effect on 
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Fig. 10. Metal temperature after three hours for 500 Hz (left) and 2650 Hz (right). 
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Fig. 11. Joule effect after five minutes for 500 Hz (left) and 2650 Hz (right). 


the crucible wall. On the contrary, when working with high frequency the power distribution 
in the furnace remains almost unaffected in the presence of molten material. 

We also present in Fig. 13 the velocity field for both frequencies. When working with low 
frequency the depth of penetration is higher, and Lorentz’s force becomes stronger than 
buoyancy forces. At high frequency, instead, the low skin depth makes Lorentz’s force almost 
negligible and buoyancy forces become dominant. This can be seen in this figure: at high 
frequency the molten metal is moving by natural convection, thus it tends to go up near 
the hot crucible, except in the upper part, probably due to the boundary condition we are 
imposing. At low frequency the electromagnetic stirring enforces the metal to go down close 
to the crucible, and a new eddy comes up in the bottom of the furnace. 

Finally, in Fig. 14 we represent the variation in time of the Joule effect in each material, along 
with the heat losses through the refrigeration tubes. In order to attain the desired power, 
higher intensities are needed working at low frequency, which causes stronger ohmic losses in 
the copper coil. Moreover, at low frequency the load begins to melt after 60 minutes, affecting 
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Fig. 12. Joule effect after three hours for 500 Hz (left) and 2650 Hz (right). 
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Fig. 13. Velocity fields after three hours for 500 Hz (left) and 2650 Hz (right). 
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Fig. 14. Joule effect and heat losses through the tubes, for 500 Hz (left) and 2650 Hz (right). 


the performance of the furnace and enforcing to increase the intensity, consequently increasing 
the power losses in the coil. When working at high frequency the performance of the furnace 
is almost unaffected by the presence of molten material. It is also remarkable that at the first 
time steps the power losses in the coil match the heat losses through the water tubes. When 
time increases the heat losses through the tubes become higher due to the heat conduction 
from the crucible across the refractory layer. 
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